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Abstract 

1 1 ^ We discuss continuous cascade models and their potential for modelling 

the energy dissipation in a turbulent flow. Continuous cascade processes, 
expressed in terms of stochastic integrals with respect to Levy bases, are 
examples of ambit processes. These models are known to reproduce exper- 
ts imentally observed properties of turbulence: The scaling and self-scaling 
C of the correlators of the energy dissipation and of the moments of the 

coarse-grained energy dissipation. We compare three models: a normal 
model, a normal inverse Gaussian model and a stable model. We show 
£h that the normal inverse Gaussian model is superior to both, the normal 

and the stable model, in terms of reproducing the distribution of the en- 
i ^ | ergy dissipation; and that the normal inverse Gaussian model is superior 

to the normal model and competitive with the stable model in terms of 
^SJ reproducing the self-scaling exponents. Furthermore, we show that the 

^ presented analysis is parsimonious in the sense that the self-scaling ex- 

ponents are predicted from the one-point distribution of the energy dis- 
sipation, and that the shape of these distributions is independent of the 
Reynolds number. 
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Since the pioneering work of Kolmogorov |30| and Oboukhov 33 , where the 
turbulent energy dissipation is assumed to be log-normally distributed, the 
small-scale intermittency of the energy dissipation in turbulence has received 



much attention 24 44 . The small scale intermittency is primarily expressed in 



terms of multifractal and universal scaling of inertial range statistics, including 



extended self-similarity 10 , scaling and self-scaling of correlators 39 , and the 



statistics of breakdown coefficients 20 



Early attemps to model the rapid variation of the turbulent velocity field 
include [TJ[9j[l2j[23j[28j[45] ( among many others). Such phenomenological ap- 
proaches are sometimes called "synthetic turbulence" and can be divided into 
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two classes. The first direction starts from modelling the velocity field and 
derives the model for the energy dissipation by taking squared small scale incre- 
ments. The second line of investigation focuses on modelling the energy dissipa- 
tion field and derives the velocity field by various, partly ad hoc, manipulations. 
The approach presented here focuses on modelling the energy dissipation as the 
fundamental field, not derived from an a priori velocity field. 

In [45| , an iterative, geometric multi-affine model for the one-dimensional 
velocity process is contructed and some of the basic, global statistical quanti- 
ties of the energy dissipation field are derived. However, this discrete, dyadic 
approach does not allow to give explicit expressions for more specific statistical 
quantities. In the approach discussed here, all main statistical quantities, like 
rt-point correlations and probability densities can be derived analytically since 
our approach defines the energy dissipation field as an explicit closed expression 
that is mathematically tractable and does not involve an iterative procedure. 

Another dyadic, iterative approach for the construction of the velocity field 
is discussed in [9j. Their model is based on a wavelet decomposition of the 
velocity field combined with a multiplicative cascading structure for the wavelet 
coefficients. The energy dissipation field is then derived from small scale in- 
crements, again leaving only limited room for an exact analytical treatment of 



higher order statistics. However, as discussed in 23 , such wavelet approaches 
are superior over discrete geometric approaches as they allow to model station- 
arity in a mathematical more rigorous way. The approach discussed here does 
not suffer from problems related to mathematical rigour and no iterative limit 
arguments are needed for the construction. A related and interesting wavelet- 



based approach is discussed in 12 , which allows for a sequential construction of 
the field. A further wavelet-based approach [l] builds on random functions and 
their orthogonal wavelet transform. The authors show that to each such random 
function there is an associated cascade on a dyadic tree of wavelet coefficients. 
The performance of the model is illustrated by numerical examples, with little 
analytical insight. 

The models 9|45 fail to incorporate skewness for the velocity increments 28 , 
a basic property of turbulent fields. As an alternative approach, [28] proposes 
a combination of a multiplicative cascade for the energy dissipation, the use of 



Kolmogorov's refined similarity hypothesis 30 and an appropriate summation 
rule for the increments to construct the velocity field. Here, again, only discrete 
iterative procedures are employed which make analytical statistical statements 
very difficult. 

It is important to mention that the continuous cascade fields discussed in 
this paper are potentially useful in the above cited works [TJ[9j[T2j[23j|28] as 
a candidate for a closed continuous and mathematically tractable and flexible 
version of the discrete multiplicative procedures employed there. 

Discrete and continuous random cascade processes have proved useful in 
describing phenomenologically the small-scale behaviour of the turbulent energy 



dissipation 11 18.25-27 31 32 38 . In 17 the surrogate energy dissipation is 



modelled as a discrete random multiplicative cascade process. Choosing the law 
of the cascade generators to be log-normal yields the Kolmogorov-Oboukhov 
model. A continuous analogue to the discrete multiplicative cascade processes 
is formulated in terms of integrals with respect to Levy bases and has been 
shown [4| [39[|4l] to be computationally tractable and to accurately describe the 
two- and three-point statistics of the energy dissipation. 
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In the cited works, focus is on the modelling of n-point statistics of the energy 



dissipation, not the distribution of the energy dissipation itself. Indeed, 39 
concludes with a remark that the law of the Levy basis driving the cascade 
model should be inferred and its dependency on the Reynolds number should 
be investigated. 

Both, discrete and continuous multiplicative cascade processes, suggest that 
the law of the logarithm of the energy dissipation should be infinitely divisi- 
ble. Infinite divisibility is necessary for the cascade models discussed here to 
be defined in a mathematical rigorous way. Furthermore, it greatly simplifies 
analytic calculations of some of the statistical properties of the models. Among 
the infinitely divisible distributions are the normal, stable, and normal inverse 
Gaussian distributions. These three classes of distributions each have their own 
tail behaviour. 

The use of stable Levy bases for modelling of the energy dissipation has been 
investigated in 20 and it is concluded by analysing the breakdown coefficients 
that "except for the log-normal limit, this leaves no room for the log-stable 
modelling of the turbulent energy cascade." The present paper investigates 
the alternative of using a normal inverse Gaussian Levy basis to model the 
energy dissipation and addresses the one-point distributions, multifractality of 
the coarse-grained energy dissipation, and two-point statistics. 

The use of normal inverse Gaussian distributions in turbulence modelling 



is not new. In 17 , the parameters of the normal inverse Gaussian distribu- 
tion cascade generator are estimated from scaling exponents and cumulants, 
which are moment estimates, notorious for their sensitivity to outliers. Indeed, 
estimation of the normal inverse Gaussian parameters (and those of other dis- 
tributions) may not be feasible from sample moments. In this paper we will 
apply likelihood methods which do not suffer from the same problems as the 
moment based methods. 

Our motivation for the use of the normal inverse Gaussian distribution is not 
based on physical arguments. We mainly exploit that these distributions are 
flexible, yet simple, and capable of attaining a wide range of shapes. Some short- 
comings of the normal inverse Gaussian distribution will be discussed. Being a 
normal mean- variance mixture with an inverse Gaussian as mixing distribution, 
the tail behavior of the mixture is related to that of the mixing distribution [8] . 
Hence one obtains a recipe for the construction of distributions with a pre- 
scribed tail behaviour. The infinite divisibility is necessary for the calculus of 
the stochastic integrals used in this paper. 



The She-Leveque-Dubrulle model 21 42 43 , which has been shown to ac 



curately predict the scaling exponents of the coarse-grained energy dissipation 
and the structure functions of the velocity increments, prescibes a Poisson dis- 
tribution for the log-energy dissipation. We find that the Poisson distribution 
provides a poor fit and that the normal inverse Gaussian distribution is a clear 
improvement. It should be noted that, being infinitely divisible, the Levy-Ito 
decomposition expresses the normal inverse Gaussian distribution as a linear 
combination of Poisson distributions. 

The paper is organised as follows. Section [2] provides some background on 
the data analysed in this paper. Section [3] recalls the construction of continuous 
cascade processes in terms of integrals with respect to Levy bases. Section [4] 
applies the theory to the data and shows how the distribution of the surro- 
gate energy dissipation determines the scaling and self-scaling exponents of the 
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Data sets no. 1 and 13 
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Figure 1: The spectral density function for data set no. 1 (lower) and no. 13 
(upper). The slope of the lines is —5/3. 

two-point correlators of the energy dissipation and the coarse-grained energy 
dissipation. Section [5] concludes. The two appendices provide necessary back- 
ground on the normal inverse Gaussian distribution and integration with respect 
to Levy bases. 



2 Background on the data 

We analyse thirteen data sets, each consisting of approximately sixteen million 
one-point time records of the velocity component in the mean stream direction in 



helium gas jet flow 15 . The time series can be assumed to be stationary. In 19 



eleven of the thirteen data sets are used in an analysis of the intermittency 



exponent of the turbulent energy cascade. In particular, 19 table I] summarises 



useful information about the data sets. Data sets no. 6 and 9 are not considered 
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Let u denote the velocity component in the mean stream direction, and let U 
denote the mean stream velocity. Since the flow can be assumed to be homoge- 
neous and isotropic, we use the surrogate energy dissipation e(x) — 15v(du/dx) 2 
as a proxy for the energy dissipation and henceforth omit the "surrogate" pred- 
icate. Here x denotes the position along the mean stream direction, and v 
denotes the viscosity. We apply the Taylor frozen flow hypothesis to express 
the energy dissipation in terms of the measured time series. We do not invoke 
an "instantaneous Taylor correction" |15| since it introduces spurious effects in 
correlators of the energy dissipation. Since any change of the energy dissipation 
by a multiplicative constant is inconsequential for the conclusions, we scale the 
energy dissipation to have unit mean. Finally, the derivative du/dt is calculated 
from the discrete samples using interpolation with third-order splines. 

The resolution scale Ace (mean velocity U times sampling time At) is 1-5 
times the Kolmogorov length -q = (v 3 / (e)) 1 / 4 , and the Taylor-microscale based 
Reynolds number ReA varies from 85 to 1181. Table [l] lists Rca in addition to 
parameters that will be explained later. 

The spectral density of the velocity component in the mean stream direc- 
tion is estimated using Welsh's overlapping segment averages with a Hanning 
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Data set no. 1 




Figure 2: Exerpts of the time series for data sets no. 1 and 13. The sampling 
resolution of data set no. 1 is approximately four times that of data set no. 13. 
The length of both exerpts is approximately 500r/. 
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taper, a block length of 500 000, and an overlap of 50%, see e.g. 34 for de- 
tails about the method. Figure [l] shows the spectral density for data sets no. 1 
and 13. The inertial ranges are found at frequencies from approximately 10 Hz 
to 100 Hz and 200 Hz to 60 kHz, respectively, after which the transition to the 
dissipation range occurs. The slope in the inertial range is approximately —5/3, 
obeying Kolmogorov's 5/3-law. For most of the data sets, the spectral density 
attains an almost constant value at the highest frequencies. We interpret this 
as instrument noise. As shown by fig. [2] (top) this does not appear to distort 
the data significantly. For some of the data sets, the resolution scale is almost 
five times the Kolmogorov length. The effect of this becomes apparent in fig. [2] 
(bottom) where the time series appears much less smooth than in fig. [2] (top) 
where the resolution scale is close to the Kolmogorov length. If data set no. 1 is 
subsampled by retaining every fourth record, a picture similar to fig. [2] (bottom) 
is obtained. Having no means to improve the resolution scale of the measured 
data sets, we consider downsampling the data to improve the smoothness so that 
calculation of the derivative du/dt becomes reliable. As will be remarked later, 
the influence of the downsampling on the shape of the one-point distribution of 
the energy dissipation is negligible and it significanly decreases the scatter of all 
quantities derived from moments estimates. Henceforth all graphs and tables 
refer to data that has been downsampled by a factor of two. We found little or 
no change when downsampling by a factor of three or four compared to a factor 
of two. 

To assess the existence of moments, we consider, for a given time series y — 
{yi, . . . , Vn}, the relative absolute sample moment M%(y)/Mff(Y) of order p 
and sample size n < N, where M%(y) = ^Y^i=i\Vi\ P - Figure ^ shows the 
relative absolute sample moments of du/dt of orders 6 and 8 for data set no. 7. 
We see that the sample moments are corrupted by outliers and that removal of 
outliers increases the reliability of the moment estimates. We will assume that 
du/dt has moments of order at least 8. The removal of outliers does not distort 
the parameters of the one-point distribution of the energy dissipation. 



3 A Levy based continuous multiplicative cas- 
cade model 

In this section we present the model for the energy dissipation. The model 
provides a link between the distribution and the two-point correlators of the 
energy dissipation. 



In 39 40 the energy dissipation e is modelled as a (1 + l)-dimensional 
stochastic process (one dimension in space and one in time) given as the expo- 
nential of an integral with respect to a Levy basis Z on R 2 , 

e(x, t) =exp( Z(dx' dt')) = exp(Z(A(x, *))), 

where A(x, t) C R 2 is called the ambit set. We will assume that the Levy basis is 
homogeneous (see appendix |B|, and that the ambit sets are defined by A(x, t) = 
(x,t) + A for some bounded set A C R 2 . Thus we ensure that e is stationary 
in space and time. The model ([I]) is an example of a random multiplicative 
cascade process in continuous space and time. Since only one-point time series 
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Data set no. 13 
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Figure 3: Relative absolute sample moment of order 6 (top) and 8 (bottom) of 
du/dt before (black) and after (gray) removal of outliers for data set no. 13. 
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are available for our analysis, we will ignore the spatial dependency and consider 
the energy dissipation as a function of time alone and at a fixed position in space, 

e(t) = cxpf / Z(dx'dt')) = exp(Z(A(t))), (1) 
VA(t) ' 

where A(t) = A + (0,t). However, it must be emphasised that the model is 
equally capable of modelling specific spatial statistics in one or more spatial 
dimensions. 

The details of the derivation of various properties of ([!]) are given in 39 40 



For convenience we repeat in the following three subsections the formulas we 
need for the subsequent analysis. 

3.1 One-point distributions 



The distribution of loge is given by (20) or (21 1. Since the integrand in ([I]) is 
constant, we have that 

K(stloge(t))=K(s$Z')vol(A), (2) 

where K(s$X) = log(exp(sJQ) denotes the logarithm of the Laplace transform 
of the random variable X, Z' denotes the Levy seed corresponding to the Levy 
basis Z, and vo\(A) denotes the volume (or area) of the set A. By pj), the 
distribution of the energy dissipation is closely related to the distribution of 
the Levy seed Z' . For brevity we define K(s) = K(s \ Z') and K(p,q) — 
K(p + q)-K(p)-K{q). 

3.2 Coarse-grained energy dissipation 

Let ei denote the coarse-grained energy dissipation, 

Y ft+l/2 

ei(t) = j e(s)ds. 

Kolmogorov's refined similarity hypothesis states that 

u(t + !)-«(*) l = V{hi) 1/3 , 

where V is a random variable independent of the energy dissipation. This 
hypothesis provides a relation between the scaling exponents C p of the moments 
((u(t + I) — u(t)) p ) ~ K p of longitudinal velocity differences and the scaling 
exponents cr p of the moments (ef ) ~ l ap of the coarse-grained energy dissipation, 
namely 

Cp =P/3 + cr p/3 . 

The exponent cr 2 is the intcrmittency exponent which quantifies the deviation 
from Kolmogorov's 1941 theory where Cp = p/3. 

From the scaling property (ef) ~ l ap we derive the self-scaling property 
(el) ~ (ef) a i/' 7 p where the ratio a q /a p is the self-scaling exponent. It is shown 
in 40 that ([T]) implies that 



C A (K(p)-pK(l)) 
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for integral values of p where Ca is a constant depending only on the size and 
shape of the ambit set A. It follows that the self-scaling exponent a q /a p is 
independent of the ambit set and hence only depends on the distribution of the 
Levy seed Z' . 



3.3 Two-point correlators 

The two-point correlator c Pj9 of order (p, q) of the energy dissipation is defined 
as 



(e(OMQ) 9 ) 

e(l)*>M<W 



The sum p + q is the total order of the correlator c Pj9 . It is shown in 39 40 
that 



Cp, q (l) = exp(K(p, q) vo\{A{l) n A(0))). 



(3) 



We observe that the exponent in ^ is expressed as a product where the first 
factor depends only on the Levy basis and the order of the correlator, and the 
second factor depends only on the overlap of the ambit sets. This provides a 
way of modelling a wide range of correlators, since the shape of the ambit set, 



under suitable assumptions, can be determined from the correlator, see 39 for 



details or subscc. [3^4] for a particular example. 
Experiments 19 41 reveal that, at least for p - 



lators exhibit scaling, 



c M (i)«r (M) , 



q < 3, the two-point corre- 



(4) 



in a range of / comparable to the inertial range of the velocity structure func- 
tions. It is straightforward to show that the scaling exponent t(1, 1) is equal 
to the intermittency exponent o-i- This follows from the fact that the second 
derivative of (l 2 ef(t)) with respect to I behaves as c\^(l) 
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In [39] it is shown that the two-point correlators also enjoy the property of 
self-scaling, 



,(0 



for an even wider range of I, in analogy to extended self-similarity 10 of struc 



ture functions. Here r(pi,qi;p2,q2) = T {P2,Q2)/T(pi,qi) is the self-scaling ex- 
ponent. 

By ([3]) we immediately obtain self-scaling of the correlators under the model ([!]) , 



,(0 



l(0 



K(pi,qi;p2,q2) 



where 



K(pi,qi;p 2 ,q2) = K(p 2 ,q2)/K(p 1 ,q 1 ) 

is the self-scaling exponent. As noted in |39j, the self-scaling property is in- 
dependent of the shape of the ambit set and thus scaling of the correlators is 
not necessary for self-scaling of the correlators. The Levy seed Z' determines 
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Figure 4: The boundary of the ambit set ^ for the values 6 = 1 (solid) ,9 = 2 
(long dashes), and 9 = 100 (short dashes). The other parameters arc T = 1 and 
L = 10. 



through ([2]) the distribution of the energy dissipation and therefore the self- 
scaling exponents. Both, the distribution and the self-scaling exponents, can be 
estimated from data and hence compared with the model. 

The model ([I]) is also able to reproduce empirical three-point correlators 41 , 
but these are not considered in the present paper. 

The correlators are moment estimates. Therefore they may not exist be- 
yond a certain order. When they exist, they are sensitive to noise and outliers, 
particularly at high orders. By Holder's inequality we have 

(e(t 1 ) p e(t 2 ) q ) < (e(ti) p+9 ) p/(p+9) (e(t 2 ) p+9 ) 9/(?,+9) . 

It follows that the correlator c p q exists provided the velocity derivative du/dt 
has finite moment of order 2(p + q). Under the model it follows by ([2| that 
c Pi9 exists if and only if the Levy seed Z' has exponential moments of order 



3.4 Ambit sets and scaling of correlators 

While the self-scaling of the correlators and the self-scaling of the coarse-grained 
energy dissipation is independent of the shape of the ambit set A, the scaling 
property Q requires A to be specified appropriately according to In this 
subsection we consider one such specification. We define for T > 0, L > 1, and 
9 > the function g by 

»(*)=( 0<«<T. (5) 



.l + (t/(T/£))« 
In the limit 9 = 00 we have that 



g(t) 



Jl 0<t<T/L, 
\T/(Lt) T/L < t <T. 



The ambit set A is given as 

A = {(x,t) \0<t<T,\x\ <g(t)}. (6) 
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Figure [4] shows examples of the ambit set for three values of the parameter 9. 
By (§, 

(l) = exp (K(p,q) J 2g(s)ds). (7) 



It follows that 

rf log Cp,g(Z) 

dlogl 



-K(p, q)2T/L, T/L < / < T. (8) 



Hence if (p, q)2T/L is the scaling exponent of c Pi9 . In the limit = oo we have 
perfect scaling, 

dl J C ™ {l) = -K{p, q)2T/L, T/L<1< T 
a log I 

The parameter 9 is merely a tuning parameters to account for imperfect scaling 
of the correlators. In view of ([3]), where the overlap of the ambit sets determines 
the correlator, T may be interpreted as the decorrelation time, and T/L is a 
time scale near and below which the scaling behaviour has terminated. 

3.5 Three Levy bases 

We consider three distributions of the Levy seed Z'\ 

{N(/i,5) normal (N), 

S(a,ft,n,5) stable (S), 
NIG(a,/?, fx, 5) normal inv. Gaussian (NIG). 

The normal Levy seed yields a log-normal model for the energy dissipation, 



consistent with the Kolmogorov-Oboukhov theory 33 . The stable Levy seed is 
included to enable comparison with |20| . The use of a normal inverse Gaussian 
Levy seed can be motivated in several ways, though none of them is based on 
physical arguments. The normal inverse Gaussian distribution is computation- 



ally tractable and allows a wide range of distributions (see fig. 11 ). Furthermore, 



the representation ( 19 1 shows that the normal inverse Gaussian distributions 
form a natural generalisation of the normal distributions by allowing the vari- 
ance of the normal distribution to be random. Since the variance is positive, 
and if there is a "typical variance" , it is natural to model the variance using an 
unimodal distribution on the positive real line. A flexible distribution of this 



type is the inverse Gaussian law as shown in fig. 12 

While other infinitely distributions could be considered as well, the three 
classes of distributions above exhaust a wide range of unimodal distributions 
with a specific tail behaviour: The normal distributions have "light tails", the 
stable distributions have "heavy tails", decaying algebraically, and the normal 
inverse Gaussian distributions have "semi-heavy tails", decaying exponentially, 



see (15). 

The common symbols for the parameters (a, ft, fi, S) have been chosen due to 
the similarity of their interpretations. In all three cases, fj, is a location param- 
eter, S is a scale parameter, a and ft are shape parameters, and ft determines 
the asymmetry of the distribution. The domain of the parameters is /i £ R and 
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5 > in all three cases; < a < 2 and —1 < f3 < 1 in the stable case; and 
\/3\ < a in the normal inverse Gaussian case. The parametrisation of the stable 



distribution is chosen to follow 37 eq. (1.1.6)], so that the log-characteristic 
function is given by 

log(exp(isS(a, /?,//, (5))) 

_ jifxs - S a \s\ a (l - i/3sign(s)tan(7ra/2)), a ^ 1, 
[i/is-£|s|(l-i/3f sign(s)log|s|), a = 1. 



It follows from (20 1 or (21 ) that the distribution of loge in each of the three 
cases is 

{N(vol(A)/i, vo\(A)S) Z 1 ~ N, 

S(a,^,vol(A)/i,vol(A) 1 /«J) Z' ~ S, (9) 
NIG(a, /?, vol(A)/i, vol(A)<5) Z' - NIG. 

We see that the shape parameters of the distribution of log e are, in all three 
cases, identical to the shape parameters of the distribution of the Levy seed Z' ', 
and that the location and scale parameters of log e are multiplied with a factor 
determined by the size of the ambit set. 

The normal distribution has exponential moments of all orders. Provided 
j3 = —1, the stable distribution has exponential moments of all non-negative 
orders. We will therefore assume that j3 = — 1 in the stable case. The normal 
inverse Gaussian distribution has exponential moments of order s if 

s<a-(3. (10) 

In these cases the log-Laplace transform of the Levy seed is given by 

Lis+\5s 2 Z'~N, 
K(s) = \ [MS- (V/cos(7ra/2))s Q S' ~ S, (11) 

[MS + 5(^a 2 - /3 2 - y/a 2 - (/3 + s) 2 ) Z' ~ NIG. 



It follows that 



'i<5 Z'~N, 

if(p,g) = /s(p,?)^-((J '/coB(7ra/2)) ^' ~ S, (12) 
5 Z' ~ NIG, 



where 



«(P, 9) = 



'pg Z' ~ N, 

(p + g) Q -p Q -q Q Z'~S, 

- ^a 2 - /3 2 - yja 2 -(P + p + q) 2 



Z' ~ NIG. 



Note that under the normal inverse Gaussian model, we must require that p+q < 
a — j3 to ensure the existence of the correlator of total order p + q. Note also 
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Data set no. 13 




Figure 5: Distribution of the logarithm of the energy dissipation represented 
by the logarithm of the probability density function for data set no. 13 down- 
sampled by a factor two: data (circles), normal fit (long dashes), stable fit 
(short dashes), normal inverse Gaussian fit (dots and dashes), and normal in- 
verse Gaussian fit constrained to possess finite exponential moments of order 
4.5 (black curve). 



that k depends only on the order (p, q) and the shape parameters. The self- 
scaling exponents of the correlators are now given by 

In a similar manner it follows that the self-scaling exponents a q /a p of the coarse- 
grained energy dissipation depend only on the order {p, q) and the shape pa- 
rameters of the Levy seed. We therefore conclude that both kinds of self-scaling 
exponents are uniquely determined by the orders under consideration and the 
shape of the one-point distribution of the energy dissipation. 



4 Data analysis 

In this section we apply the model from section|3]to the thirteen data sets. The 
one-point distribution of the energy dissipation is used to predict the scaling 
and self-scaling exponents of the two-point correlators and the coarse-grained 
energy dissipation. 



4.1 One-point distributions 

Figure [5] shows that the distribution of log e has a distinct, non-Gaussian, asym- 
metric shape. The shape is independent of the downsamling applied to smooth 
the data. Four parametric distributions are fitted to the data using likelihood 
methods: normal, stable, normal inverse Gaussian, and normal inverse Gaussian 
constrained to have finite exponential moment of order 4.5 according to (10 1. 
Clearly, the normal inverse Gaussian distribution provides an excellent fit for 
all amplitudes. The left tail is overestimated by the stable distribution. The 
unconstrained normal inverse Gaussian distribution is also able to accurately 
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1.40 


4.49 


0.129 


0.130 


8 


463 


2.51 


-1.99 


2.54 


3.31 


1 


67 


-1.92 


1.39 


3.76 


0.128 


0.092 


9 


703 


2.50 


-2.00 


2.38 


3.03 


1 


65 


-1.86 


1.34 


4.91 


0.108 




10 


885 


2.50 


-2.00 


2.45 


3.13 


1 


65 


-1.89 


1.37 


6.72 


0.107 


0.089 


11 


929 


2.50 


-2.00 


2.44 


3.11 


1 


65 


-1.88 


1.36 


7.16 


0.112 


0.079 


12 


985 


2.50 


-2.00 


2.43 


3.07 


1 


65 


-1.89 


1.36 


9.53 


0.104 


0.105 


13 


1181 


2.49 


-2.01 


2.30 


2.88 


1 


63 


-1.83 


1.33 


9.72 


0.097 


0.061 



Table 1: Summary of the thirteen data sets. The Taylor micro-scale Reynolds 
numbers are from [15] . The normal inverse Gaussian parameters (constrained 
to yield finite exponential moment of order 4.5 according to (10)) and the sta- 
ble parameters are estimated from the distribution (|9| of log e using likelihood 
methods. Aa; denotes in units of the Kolmogorov length rj the effective spatial 
resolution after downsampling by a factor of two. The effect of downsampling by 
a factor of two changes the normal inverse Gaussian shape parameters a and f3 
by around 0.4 % and the stable shape parameter a by around 2%. r(l, 1) is the 
intermittency exponent estimated by fitting ^ and ^ to the correlator c^i, 
and t'(1, 1) is the intermittency exponent determined in [19|. 



fit the distribution, though it slightly underestimates the left tail. In what fol- 
lows, any reference to the normal inverse Gaussian distribution will imply the 
constrained version. 

It is necessary to constrain the shape parameters of the normal inverse Gaus- 
sian distribution to ensure the existence of the correlators up to a certain total 
order. This constrainment is a delicate issue. Firstly, it introduces the problem 
of determining an appropriate upper order. And secondly, it implies a slightly 
poorer fit. The choice of 4.5 as the upper order is somewhat arbitrary but pro- 
vides a compromise between fitting the distribution of the energy dissipation 
accurately and predicting the self-scaling exponents of the correlators and the 
coarse-grained energy dissipation. We conjecture the existence of a more ap- 
propropriate infinitely divisible distribution, perhaps within the class of normal 
mean- variance mixtures. For the present paper, the normal inverse Gaussian 
will suffice. 

Table [T] summarises the estimated parameters for each of the thirteen data 
sets. Note that for both, the normal inverse Gaussian distribution and the sta- 
ble distribution, the estimated shape parameters do not depend on the Taylor 
micro-scale Reynolds number, while the location and scale parameters vary only 
slightly. This is a clear indication of universality of the distribution of the energy 
dissipation, at least within the thirteen data sets considered here. Downsam- 
pling the data by a factor of two changes the normal inverse Gaussian shape 
parameters a and /3 by only around 0.4% and the stable shape parameter a by 
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Data set no. 7 




1.5 2.0 3.0 



Figure 6: Scaling and self-scaling of the moments (ef) of the coarse-grained 
energy dissipation for data set no. 7. Top: (e^) as a function of the lag I. The 
scaling behaviour is far from pronounced. Bottom: (e?) for p — 3, 4 as a function 
of (ef). The self-scaling behaviour is clear. 

around 2%. 

4.2 Coarse-grained energy dissipation 

The behaviour of the moments (ef ) as a function of the lag I does not reveal 
any pronounced scaling range, see fig. [6] (top) as an example. Compare with 
fig. [8] which demonstrates a clear scaling range for the two-point correlators. 
However, as shown by fig. [6] (bottom) the moments display a clear self-scaling 
behaviour, and the self-scaling exponents a q /a p can be extracted. 

The estimated self-scaling exponents can now be compared with the predic- 
tions of the model ([l]) and with the She-Leveque-Dubrulle model. The former 
predicts 

= K(q) - qK(l) 
u v K(p)-pK{lY 

while the latter predicts 

<jq _ -2 g /3 + 2(l + (2/3)9) 
<j p ~ -2p/3 + 2(1 + (2/3)p) ' 
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Figure 7: Self-scaling exponents (73/0-2 (top) and 0-4/0-2 (bottom) of the moments 
(ef) of the coarse-grained energy dissipation: data (circles), normal model (long 
dashes), stable model (short dashes), normal inverse Gaussian model (solid), 
and She-Leveque-Dubrulle model (dots and dashes). 
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Figure 8: The two-point correlator Cn as a function of lag I for data set no. 7 
in double logarithmic representation (circles). The solid curve shows the least 
squares fit of ^ in the double logarithmic representation. The estimated scaling 
exponent r(l, 1), decorrelation time T, and auxiliary parameter 8 are shown in 
the figure. 



Figure [7] (top) shows that the normal inverse Gaussian and the stable model 
agree well with both data and the She-Leveque-Dubrulle model for the self- 
scaling exponent 0-3/172 ■ For 04/02 fig. [7] (bottom) displays a somewhat greater 
discrepancy between the models. The She-Leveque-Dubrulle model agrees best 
with data for the data sets with the lowest Re> while the stable model agrees best 
for the other data sets. The normal inverse Gaussian model slightly overpredicts 
the self-scaling exponents, but in general it provides a prediction of comparable 
accuracy. 

The normal inverse Gaussian model can be forced to align better with the 
stable model and the She-Leveque-Dubrulle model by requiring exponential mo- 
ments of order 5 to exist. This, however, compromises the prediction of the 
self-scaling exponents of the two-point correlators (subsec. 4.4). 



4.3 Two-point correlators: scaling 

Since the sample moments of the velocity derivative become increasingly unre- 
liable with increasing order, we consider only the orders (p, q) for which p and q 
are positive half-integers with p + q < 3. This leads to 36 non-trivial combina- 
tions in the analysis of self-scaling of the correlators. The self-scaling exponents 
satisfy 

r(pi,qv,P3,q3) = r(pi, qi;p 2 , g , 2)r(p 2 , <72!P3, 93), 

and it is therefore sufficient to consider eight combinations. 

Figure [8] shows that the correlator ci t i for data set no. 7 exhibits scaling for 
lOAi < I < 200Ai. The figure also shows that the correlator Q determined 
by ([5]) provides a very good fit to the data for lags / > lOAt. The parameter 
9 determines the behaviour of the correlator for lags close to the decorrelation 
time T. At lags I < 10 At, the correlators are corrupted by surrogacy effects [i~6] . 
Lags / < lOAt have therefore been excluded from the fit. Similar results hold 
for the other data sets and other orders of the correlators. 
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The parameter L, which determines the extend of the scaling range of the 
correlators, cannot be identified since surrogacy effects corrupt the correlator 
at small lags. We have chosen L = 10 4 as it allows for clear scaling at all lags 
lAi < I < TAt for all thirteen data sets. 

The scaling exponents, under model ([I]), follow from ([8]), (J9j) , and (12 1, 



l<5 e Z' ~ N, 

— s a 

r(P^) = CL,e-n( P , q )-{ cos ^ a/2) Z'~S, (13) 
J e Z'-NIG, 

where Cl,b is a factor which depends only on L and 9. Since L cannot be reliably 



identified by the data available, ( 13 ) does not allow us to predict the correlator 
scaling exponents from the distribution of the energy dissipation alone. 

The scaling exponents may still be estimated directly by fitting §5§ and ^ 
to the correlators. Table [I] shows the estimated intermittency exponent t(1, 1) 
where it is also compared to the estimate found in [19] . The two estimates are 
in reasonable agreement. The differences are likely due to the differences in the 
estimation procedures. We note that for all the data sets, the intermittency ex- 
ponent is approximately half of the prediction under the She-Leveque-Dubrullc 
model where r(l, 1) = 2/9. 

As a final consistency check, we note that, according to ([8]), the scaling 
exponents satisfy a relation of the form 

r(p, q) = a p+q -a p - a q . (14) 



By (11 1, we may choose a\ = 0. The other unknowns are found by solving the 



linear least squares system corresponding to (14). The root mean square of the 



relative error (a p+q — a p — a q ) /r(p, q) — 1 is below 1 % except for data sets no. 



1, 2, and 5 where it is 3%, 9%, and 8%, respectively. We conclude that (14) is 
fulfilled to a satisfying degree. 



4.4 Two-point correlators: self-scaling 

The self-scaling of the correlators is confirmed by fig. [9] Under model (JlJ , the 
self-scaling exponent r(pi, <72;f>2, (tit) is given by the ratio k(p2,<72)/k(pi,<Zi)- It 
depends only on the shape parameters of the distribution of log e and on the 
order of the involved correlators. In particular, the self-scaling exponents are 
independent of the ambit set and the location and scale parameters of log e. 
Therefore, the self-scaling exponents are predicted directly from the estimated 
shape parameters listed in tab. [I] and they inherit the ReA-independence from 
loge. 

Figure [T0| shows the estimated self-scaling exponents for all the data sets. In 
most cases the predictions of the normal inverse Gaussian model and the stable 
model arc within the variation of the data point, though the normal inverse 
Gaussian model tends to predict larger values than the stable model. Both 
models are clearly superior to the normal model. The stable model appears 
to be slightly more accurate than the normal inverse Gaussian model. The 
case (1.0, 1.5; 0.5, 2.5) is exceptional as all models fail to predict the observed 
self-scaling exponents by approximately 10%. 
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Figure 9: The two-point correlator c P2 ^ 2 (y-axis) as a function of the two-point 
correlator c P1iIJi (x-axis) for data set no. 7 in double logarithmic representation. 
The black line indicates the self-scaling behaviour c P2i92 = (cp lJ , ^1 ) T ( pl ' <^l;P2 ' <^2 ) 
as shown in each subfigure. Similar results hold for the other data sets. 
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Figure 10: The self-scaling exponents r(pi, qi;P2, Q2) (y-axis) for each of the 
thirteen data sets: data (circles), normal model (long dashes), stable model 
(short dashes), normal inverse Gaussian model (solid). 
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At low Reynolds numbers, the self-scaling exponents tend to deviate from the 
self-scaling exponents at higher Reynolds numbers. At low Reynolds numbers 
the inertial range is shorter. Whether this affects the estimation of the self- 
scaling exponents remains to be investigated. 

Further observations may be drawn from fig. |10| Firstly, the predicted self- 
scaling exponents are in all three cases almost independent of the Taylor micro- 
scale Reynolds number. This follows from the fact that the estimated shape 
parameters of the distribution of loge are almost independent of Re a (tab. [IJ. 
As noted previously, this independence hints at a universal character of the self- 
scaling exponents and the distribution of loge. Secondly, the shape parameters 
are derived from one-point statistics, while the self-scaling exponents are derived 
from two-point statistics. The predictability of the latter from the former hints 
at a parsimonious description of the energy dissipation as the exponential of an 
integral with respect to a Levy basis: Higher order statistics are reasonably well 
predicted from lower order statistics. 

5 Conclusion and outlook 

Thirteen time series of one-point measurements of the velocity component in 
the mean stream direction in a helium jet are analysed from the point of view 
of the (surrogate) energy dissipation. The distribution of the logarithm of the 
energy dissipation is shown to be well approximated by normal inverse Gaussian 
distributions, a property possessed by neither the normal distribution nor the 
stable distribution. Furthermore, the shape of the distribution is apparently 
universal, i.e., it does not depend on the Reynolds number, at least within the 
considered data sets. 

The two-point correlators of the energy dissipation show scaling and self- 
scaling. By modelling the energy dissipation as the exponential of an integral 
with respect to a Levy basis, a connection is established between the shape of 
the distribution of the logarithm of the energy dissipation and the scaling and 
self-scaling exponents of the two-point correlators. Thus the model is parsimo- 
nious: the two-point correlators are all, to good accuracy, determined by the 
one-point distribution. The normal inverse Gaussian distributions are compared 
to the stable distributions. While the stable distributions are also capable of 
modelling the scaling and self-scaling exponents, they are far from capable of 
accurately modelling the distribution of the energy dissipation (fig. [5| . In par- 
ticular, the stable model implies infinite variance of log e, yet the data suggest 
that loge has finite moments at least up to order 6. The use of normal inverse 
Gaussian distributions allows accurate modelling of both the correlators and the 
distribution of the energy dissipation. 

The model also allows prediction of the self-scaling exponents of the moments 
of the coarse-grained energy dissipation. Comparison to the predictions of the 
She-Leveque-Dubrulle model is made. 

While the normal inverse Gaussian distributions allow fairly accurate mod- 
elling of the one-point distribution, two-point correlators of the energy dissipa- 
tion, and the self-scaling exponents of the moments of the coarse-grained energy 
dissipation, a slight disadvantage related to exponential moments is identified. 
We conjecture the existence of a more appropriate infinitely divisible distribu- 
tion without this disadvantage. 
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The energy dissipation exhibits stylised features beyond the scaling and 
self-scaling of the two-point correlators. In [20], breakdown coefficients and 
Kramers-Moyal coefficients are employed, in particular, to evaluate the use of 
stable distributions in the modelling. A similar analysis remains to be performed 
in the case of the normal inverse Gaussian distributions. 



A The normal inverse Gaussian distribution 

The normal inverse Gaussian distributions form a four-parameter family of prob- 
ability distributions on the real line. They are a special case of the generalised 
hyperbolic distributions introduced in [5] to describe the law of the logarithm 
of the size of sand particles (see also [2[|3j|6j[T3j[22]). The generalised hyper- 
bolic distributions are applied in many areas of science, see e.g. [7 35 and the 
references therein. 

The probability density function of a normal inverse Gaussian distribution 
is given by 



where 7 = a 2 — (3 2 7 q(x) = y/1 + x 2 , and K\ denotes the modified Bessel function 
of the second kind with index 1. The real parameter p, determines the location, 
and the positive parameter <5 determines the scale. The parameters a and /3 are 
shape parameters and lie within the shape cone: \0\ < a. 

From the asymptotic property K\(x) oc x~ 1 / 2 e~ x as x —¥ 00 it follows that 
the NIG distribution has semi-heavy tails, specifically 

P df NiG(a,,3,M,<5) ( x ) ^ M~ 3/2 exp(-a|a;| + f3x) (15) 

as x —¥ ±00. This illuminates the role of a and j3 in determining the tails of 
the distribution. 

The cumulant function K(t;a,f3,p,,5) — log(exp(£A")) of a random vari- 
able X with distribution NIG (a, /?, pt, S) is given by 

K{t; a, (3, p,, S) = pt + S(j- ^a 2 -((3 + t) 2 ) , (16) 

and the radius of convergence for the cumulant function is \a — /3|. 

By differentiating ( 16 ) we obtain the following expressions for the first four 
cumulants, 



P _ 6 1 

Kl " M + d (l-p2)l/2' K2 ~ „ (1-^)3/2- 

5 p _ o 5 1 + 4p 2 



(17) 



a 2 (l_ p 2)5/2' 4 a 3 (l-p 2 ) 7 /2' 

where p = ft /a. Hence, the standardised third and fourth cumulants are 
^3 = o P *± = o 1 + 4p 2 

K 3 / 2 (fc(l-/7 2 ) 1 /2)l/2' K 2 Sail-p 2 ) 1 / 2 ' 



(18) 
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Figure 11: The normal inverse Gaussian shape triangle. The black dots mark 
positions (x, £) m the shape triangle, and the small graph nearby shows the log- 
arithm of the probability density function with the corresponding asymmetry \, 
steepness £, zero mean, and unit variance, plotted on the interval [—3,3]. 



Equations (17 1 and (18) further illuminate the roles of the four parameters: \x 
and 6 determine location and scale, respectively; j3 is related to the skewness, 
specifically the tail asymmetry (if /3 = 0, the distribution is symmetric); and a 
is related to the kurtosis. 



It follows immediately from ( 16 ) that if X\ , . . . , X n are independent normal 
inverse Gaussian variables with common parameters a and /3, but with individ- 
ual location and scale parameters fa and Si (i = 1, . . . , n), then the distribution 

of the sum X + = X± -\ H X n is NIG(a, /3, 5+) where f/,+ = + • • • + /i„ 

and 5+ = 5i + ■ ■ ■ + 5 n . Therefore, the normal inverse Gaussian distributions 
are infinitely divisible, see also [3]. 

It is often desirable to describe the NIG distributions in terms of location- 
scale invariant parameters. By letting a — da and f3 — 6/3 we have that a and /3 
are invariant under change of location and scale. For the purpose of interpreting 
the shape parameters it is sometimes advantageous to express the shape in terms 
of the location-scale invariant steepness £ and asymmetry x, defined by 

c 1 c P 



where 7 = #7 = 8\J a 1 + ft 2 , and p = (3/a = (3/a is the alternative asymmetry. 
These parameters are within the shape triangle, defined by 

{(x,oio<e<i,-e<x<a 

Figure [TT] shows the shape of the normal inverse Gaussian distributions for 
various values of the asymmetry x and steepness £. A wide range of shapes is 
possible. (See [13] for details on the shape of the family of generalised hyperbolic 
distributions) . 

If the variance is held constant, the normal distribution is obtained in the 
limit of zero steepness. For details and examples of other limiting distributions, 
see [2}[22]. 

A useful property of the normal inverse Gaussian distribution is the repre- 
sentation in terms of a mean-variance mixture of a normal distribution with 
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x 

Figure 12: The logarithm of the probability density function of the inverse 
Gaussian distribution with parameters chosen to yield normal inverse Gaussian 
distributions with zero mean, unit variance, zero asymmetry, and steepness 
£ = 0.1 (black), £ = 0.5 (long dashes), and £ = 0.9 (short dashes), see fig. 11 



the mixing distribution being an inverse Gaussian distribution, hence the name. 
Specifically, if X is normal with mean \i + fia 2 and variance a 2 , and if a 2 is 
endowed with an independent inverse Gaussian distribution with parameters 7 
and (5, then X follows a normal inverse Gaussian distribution with parameters 
(a, P, fi, 6). In short we may write 

X ~ fi + 13a 2 + all (19) 

where U is a normal distribution with zero mean and unit variance. For refer- 
ence, the probability density function of an inverse Gaussian distribution with 
parameters 7 and S is 



pdfiG( 7 ,*)(z) = ^Lx- 3 / 2 expH(x 7 2 + x~H 2 )) 



for x > 0. Figure [12] shows the logarithm of the probability density functions of 
the inverse Gaussian distribution corresponding to the three symmetric normal 
inverse Gaussian distributions in fig. |11| A s the steepness increases, the proba- 
bility that the random variance a 2 in ( |19| will attain large values increases. 

To estimate the normal inverse Gaussian parameters from data one may 
apply maximum (or pseudo) likelihood methods. The maximum likelihood esti- 
mation is non-trivial since the likelihood function is very flat near the optimum. 



The computer program "hyp" 14 implements numerical maximisation of the 
likelihood function, but, in general, non-linear optimisation algorithms may also 
be applied. The approach of expectation-maximisation developed in |29| may 
be used in conjunction with other optimisation algorithms. 



B Integration with respect to Levy bases 

The stochastic processes to be considered in the present paper are expressed in 
terms of integrals of deterministic functions with respect to Levy bases on R 2 . 
A Levy basis Z on R 2 is an infinitely divisible, independently scattered random 
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measure, i.e., to each bounded Borel subset A of R 2 , an infinitely divisible 
random variable Z(A) is associated, the random variables associated to disjoint 
subsets are independent, and the random variable associated to a disjoint union 
is almost surely equal to the sum of the random variables associated to each 
subset, provided the union is a bounded Borel subset. For more details and 



mathematical rigour we refer to 36 



The stochastic integral J f dZ of a deterministic function / : R 2 — ► R with 
respect to a Levy basis Z is defined in three steps. Firstly, for an indicator 
function 1^ we define J 1a dZ — Z(A). Secondly, by requiring that the integral 
is linear in the integrand, the integral is extended to simple functions, i.e., 
linear combinations of indicator functions. Finally, since a measurable function 
/ : R 2 — > R may be approximated by a sequence of simple functions, the integral 
J f dZ is defined to be the limit in probability of the sequence of integrals of 
the simple functions, provided this limit exists. 

An important class of Levy bases has the property that, informally, the 
distribution of the random variable associated to a subset does not depend on 
the location of the subset. In this case, we have the following fundamental 
representations of the Laplace transform and the characteristic function of the 
integral of a deterministic function / with respect to a Levy basis Z. Let 
K(s\X) = log(exp(sX)) and C(s\X) = log(exp(isX)) denote the logarithm of 
the Laplace transform and the characteristic function of the random variable X, 
respectively. Then 



K(s % j /(a) Z{da)) = J K(sf(a) J Z') da, (20) 
C(s \ [ f(a) Z{da)) = [ C(sf(a) X Z') da, (21) 



where Z' is a random variable (called the Levy seed) whose cumulant function 
is related to the Levy basis Z by 

K{s\Z(da)) =K{s%Z')da, 
C{sXZ{da)) = C{sXZ')da, 

see g for more details. It follows that the distribution of the stochastic integral 
is determined by the function / and the log-characteristic function of the Levy 
seed Z'. 
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